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I. INTRODUCTION 

Quantum Monte Carlo (QMC) simulations have be- 
come the method of choice for studying large equilibrium 
quantum many-body systems without approximations in 
more than one dimension [l], Q (in one dimension, the 
Density Matrix Rcnormalization Group has proven to be 
an extremely powerful tool 0,0]). While for small system 
sizes one may employ exact-diagonalization techniques, 
for larger ones, QMC methods provide in most cases the 
only numerical method available for exact numerical in- 
vestigation. However, employing Monte Carlo techniques 
comes at a cost. As the name itself might indicate, QMC 
methods are stochastic in nature as they are based on 
sampling the exponential number of states of the Hilbcrt 
space of the system, and there are therefore statistical 
errors associated with every measured quantity. 

QMC methods arc usually considered ideal for mea- 
surements of ground-state properties or for the determi- 
nation of thermodynamic properties of physical systems, 
as these can usually be measured to a high degree of ac- 
curacy, i.e., with very small statistical errors. While this 
is true, QMC methods are also known to be less suited 
for extracting information about excited states, which 
tend to be rather cumbersome to obtain and are typi- 
cally measured with much less accuracy. 

Excited states play a central role in many areas of 
physics and chemistry of many-body systems. Among 
these are critical phenomena and phase transitions in 
Condensed Matter Physics [f|, mass gap calculations in 
High Energy Physics various calculations pertain- 

ing to the properties of nuclei in Nuclear Physics 0, [Io| 
and the vibrational modes of large molecules in Chem- 
istry [ll], Qjj], to mention some diverse examples. Natu- 
rally, numerous attempts have been made to utilize quan- 
tum Monte Carlo methods to compute excited-state ener- 
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gies as well. However, this has turned out to be a difficult 
task because obtaining information about excited states 
involves 'isolating' specific regions within the spectrum of 
the Hamiltonian - something which can not be done by 
simple measurements of thermodynamic properties (ex- 
cept maybe the ground state at ultra- low temperatures) . 

As quantum Monte Carlo methods have evolved over 
the years, techniques to extract information about ex- 
cited states have been continuously developed. Most of 
these were based on some form of analysis performed on 
measurements of imaginary-time correlation functions as 
these provide indirect access to the spectrum of the sys- 
tem Hamiltonian (l3l - [l5| (this will be explained in more 
detail in the next section). Despite these elaborate ma- 
nipulations on the measured data, an accurate prediction 
of even the lowest excitation energies still remains a chal- 
lenge - and under some circumstances an impossibility - 
due to the large statistical errors associated with corre- 
lations with large imaginary-time differences, although 
methods to reduce these errors in certain cases through 
the use of improved estimators in cluster-based QMC 
methods [l6| or the use of 'smeared' operators in lattice 
gauge theories [I?], HU have been devised. 

In what follows, we propose a remedy to these difficul- 
ties by suggesting a way to partially optimize the manner 
in which excited-state energies, specifically the gap to the 
first excited state, are calculated from imaginary-time 
correlation functions. We do this by providing a pre- 
scription to find and then measure the most suitable cor- 
relation function available for this purpose (within some 
stated limitations). The method we suggest here is based 
on finding the operator whose imaginary-time correlation 
function is optimal for the extraction of excited-state en- 
ergies, where the optimization is based on the maximiza- 
tion of the integrated susceptibility within a space of 'ba- 
sic operators' and under appropriate constraints. As we 
shall demonstrate, this type of optimization removes, or 
at least substantially reduces, some of the difficulties as- 
sociated with dealing with the large statistical errors that 
characterize correlations with large imaginary-time sep- 
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arations. 

The paper is organized as follows. In Sec. [H] we discuss 
in some detail the basics of extracting excited-state ener- 
gies from imaginary-time correlation functions, focusing 
on the extraction of the gap to the first excited state. 
We also list some of the difficulties involved in doing so. 
In Sec. IIIII we shall present a method to find a measur- 
able operator whose imaginary-time correlation function 
is optimal for the extraction of the gap. Wc provide sev- 
eral illustrative examples of the method in Sec. IIVI and 
summarize the results in Sec. |V] along with some conclu- 
sions. 



II. ACCESSING THE EXCITATION GAP - 
CORRELATION FUNCTIONS IN IMAGINARY 
TIME 

Let us consider a many-body system described by 
the Hamiltonian H (imagine, say, an iV-body system 
of interacting spin-1/2 particles) at inverse-temperature 
j3 = 1/T (in our units, ks = 1). The thermal averages of 
physical observables are given by [3l| : 

(6) = i X Tr[6e"^] , (1) 

where O is the operator associated with the physical ob- 
servable in question. Here, Z is the partition function 
Z = Tr[c~ ,3ff ] and Tr is the trace operation. 

If the system is small enough, these thermal averages 
may be computed rather easily by exact-diagonalization 
techniques, with which the full matrix of the Hamilto- 
nian is spectrally decomposed. In this case, excited-state 
energies would simply be obtained by subtracting the 
lowest eigenvalue of the Hamiltonian matrix from other 
eigenvalues. 

For bigger systems however, where exact- 
diagonalization methods are unfeasible, one must 
almost always resort to QMC techniques to obtain 
accurate results. While the thermal averages of most 
operators of potential interest are usually very easy to 
obtain within QMC simulations (as discussed in the 
Introduction) , in order to evaluate excited-state energies 
of the system, only somewhat-indirect methods are 
available. 

The first step toward finding excited-state energies is 
the calculation of the thermal averages of different-time 
correlations of measurable operators that do not com- 
mute with the Hamiltonian (and hence are not conserved 
in time). Consider the imaginary-time 'two-point' corre- 
lation of the measurable operator O, namely (O(t)O(O)), 
where r is the imaginary-time coordinate. This expres- 
sion may be expanded in the eigen-energy basis to give: 

(O(r)O(O)) = <e^O(0)e-^6(0)) = (2) 
(E c ~^) x E \(n\d\m)\ 2 e-^- E ^ T e-^ , 

\k=0 I n,m=0 



where {E n } and {|n)} are the eigenvalues and matching 
cigenstates of the Hamiltonian H. Now, if f3 is chosen 
such that (5AEi ^> 1 where AE n = E n — Eq is the gap 
to the n-th excited state (and in particular the excitation 
gap is AEi), it is expected that the system will eventually 
relax to its ground state |0) [HI]. Under this condition, 
which wc shall assume to hold henceforth, one could de- 
fine the following correlation function with the associated 
series expansion: 

C 6 (t) = (O(r)O(O)) - (O) 2 (3) 
« ^|<0|O|n)| 2 (e- Ai ^+e- AiJ "^) . 

n = l 

Information about excited-state energies (or equiva- 
lently the gaps AE n ) is usually extracted by fitting mea- 
surement data of the correlation function or some trans- 
formation thereof, to an expression similar to the sum 
in the above equation, where usually the free parameters 
of such fits correspond to the energy gaps AE n and the 
matrix elements |(0|O|n)|. For obvious reasons, finding 
more than a few excited-state gaps is unfeasible because 
of the exponential number of free parameters involved in 
the fit and the finite number of available uncorrelatcd 
measurement data points (although attempts to model 
the some of the spectrum by a continuum of states has 
also been suggested 

Here, we shall focus the discussion on a rather basic 
type of analysis of the imaginary-time correlation func- 
tion data with which the gap to the first excited state is 
extracted, although it should be noted that more sophis- 
ticated methods of analysis exist and may be employed 
just as easily. In fact, these methods are expected in 
most cases to perform better than the simple analysis 
and provide better estimates of the excitation gap and 
possibly also limited information on higher energy levels 
[131 — 1151 ] - However, application of these methods of anal- 
ysis and comparison between them is complementary to 
the discussion here and will therefore remain outside the 
scope of this paper. 

Examination of the form of the correlation function 
given in Eq. ([3]) suggests that it might be possible to 
extract the excitation gap by analyzing the behavior of 
the correlation function at long imaginary times where 
the slowest-decaying exponent dominates the series and 
as a result the correlation function may be approximated 
in that region by: 

C7 6 (r)«|(0|6|l)| 2 e- A ^. (4) 

In this case, the simplest and most straightforward 
method of analysis for extracting the gap AEi would 
simply be fitting the logarithm of the obtained measure- 
ment data of Cq{t) acquired in the simulation with a 
straight line in the said region. For the gap to be obtained 
accurately, however, in addition to being dominated by 
the slowest-decaying term, the correlation function must 
also have small relative statistical errors, i.e., a large 
signal-to-noise ratio - which is usually a feature of short 
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imaginary-time correlations. Normally then, one looks 
for an intermediate region between r = and t = (3/2 
(the correlation function is symmetric about ft/ 2) that 
satisfies both of the above demands. It should also be 
noted that the choice of ft also plays a role in finding an 
appropriate region: If ft is chosen to be too small, there 
will be no region where only the slowest-decaying expo- 
nent survives. On the other hand, if ft is chosen to be 
too big, the system will take longer to equilibrate. 

An illustrative example of the method is given in Fig.Q] 
In the figure, results of the above analysis applied to a 
system of 64 interacting spin- 1/2 particles encoding a 64- 
bit l-in-3SAT problem augmented by transverse fields is 
presented. (The specific structure of the Hamiltonian of 
the system is discussed in Sec. IIV1 ) The figure shows the 
correlation function of the diagonal part of the Hamil- 
tonian of the system accompanied by a linear fit in an 
intermediate region aimed at evaluating the gap of the 
system. 

As the figure indicates, choosing the appropriate re- 
gion in imaginary time is not always a simple task: At 
small t, all the exponents in the sum given in Eq. ([3]) 
are expected to contribute significantly to Cq (t) in pro- 
portion to the square of their respective matrix elements 
(0|O|n}| 2 . At the other end of the range, at times closer 
to ft/ 2, it is very likely that the slowest-decaying expo- 
nent will be the only surviving one (provided that ft is 
large enough), however in practice, it may also be very 
difficult to obtain an accurate estimate of the correlation 
function there, since the signal-to-noise ratio of correla- 
tions with significant time differences is typically very 
small due to the exponential decay of the function (this 
is also evident in Fig. [J). 

Possible improvements over the above method of anal- 
ysis would involve for example, fitting the correlation 
function with hyperbolic cosines which would account 
for the signal coming from ft/2 < r < ft, or adjusting 
the fit to partially account for contributions of higher- 
energy excitations. These methods will not be discussed 
further here. 



III. OPTIMIZED CORRELATION FUNCTIONS 

In the previous section, we saw that our capability of 
accurately extracting the gap depends on the existence of 
a region where the correlation function can be nicely fit 
with a straight line (on a log-linear scale). The expres- 
sion for the correlation function given in Eq. Q shows 
that for a given value of ft, the exact shape of the cor- 
relation function is dependent upon two sets of values. 
The first is the spectrum of H which is a given property 
of the system. The second set of values is the matrix 
elements (0|O|n)| 2 which on the other hand can be par- 
tially manipulated by different choices of the operator 
O whose different-time correlations are being measured. 
The general guideline for choosing the most suitable op- 
erator to measure the correlation function of, for the ex- 
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FIG. 1: A log-linear plot of a time dependent correlation func- 
tion for an N = 64 spins system of one instance of the 1-in- 
3SAT problem with ft = 1024 (further details of the prob- 
lem can be found in Sec. IIV[I . While at early times there 
are contributions from several energy levels, there is a re- 
gion where the correlation function may be fit with a simple 
straight line (on a log-linear scale) from which the energy gap 
which is the negative of the slope could be obtained (giving 
here AEi = 0.0364). At still later times where the correla- 
tion function is significantly attenuated, the statistical errors 
become huge. 



traction of the gap would naturally be based on bringing 
the matrix element |(0|O|l)| 2 to a maximum while keep- 
ing the other matrix elements (0|O|n}| 2 for n > 1 con- 
strained. Doing so would have a two-pronged effect on 
the correlation function: it would yield a more dominant 
slowest-decaying exponent at short times and will also 
make the correlations stronger throughout. This would 
in turn substantially reduce the statistical errors at in- 
termediate and long times. 



A. Maximizing the integrated susceptibility 

The matrix elements |(0|O|n)| discussed above can 
not however be accessed or manipulated directly within 
quantum Monte Carlo simulations as this requires knowl- 
edge of the excited states of the system - knowledge that 
one does not have when the system is in its ground state. 

A question then arises as to how one could determine 
which measurable operator O has the optimal imaginary- 
time correlation function for the extraction of the excita- 
tion gap? As it turns out, there is a measurable thermo- 
dynamic physical quantity associated with every operator 
that may provide indirect access to these matrix elements 
and as wc shall sec will prove to be key in finding a well- 
suited operator for the extraction of the gap [20|. This 
quantity is the 'integrated susceptibility' of the operator 
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which is defined by: 



Xd 



E 



C 6 (r)dr 
\{0\O\n)\ 2 



(5) 



x 2(1 



-AE n /3 



The above quantity has two very favorable properties. 
Firstly, it is a zero-frequency quantity, and while it in- 
volves integration over the entire range of imaginary 
time, it can still be very easily and efficiently measured in 
the course of a simulation within a variety of QMC algo- 
rithms (see Appendix [5] for a description of its measure- 
ment within the Stochastic Series Expansion algorithm). 
Secondly, as the sum in the above equation indicates, the 
integrated susceptibility may be viewed as an estimator 
or measure of the (squared) matrix element to the first 
excited state |(0|O|l)| 2 : In particular, in cases where the 
second and higher excited states have considerably higher 
energies than that of the first excited state, Xd cou ld be 
approximated by 



Xd 



l(Q|Q|i)l 2 



(0) 



and one could therefore use the integrated susceptibility 
as an indication for the magnitude of (0|O|l)| 2 . The in- 
tegrated susceptibility may thus be used as a 'figure of 
merit' for the effectiveness of any candidate Cq(t): the 
larger Xd ' 1S > the better Cq(t) would be for extracting the 
gap. Graphically, maximizing Xq corresponds to 'lifting' 
the correlation function curve as much as possible above 
the horizontal axis thereby maximizing the area under- 
neath in the region [0,/3]. 

Given the above discussion, we are now at a point 
where we can reformulate in a more concrete manner 
the question of finding the optimal correlation function 
C 6 (t) = (O(r)6(0))-(6) 2 for the extraction of the exci- 
tation gap: Suppose that within a QMC simulation, there 
exists a set of M basic observables {Ai, . . . , Ai, . . . , Am} 
which can be easily measured in the course of the simu- 
lation [33j]. What would be the operator O = J2iLi a iAi, 
where ai are real-valued coefficients, such that \d ^ s max ~ 
imal? 

Expressing the integrated susceptibility of the operator 
O in terms of the coefficients a, , we have: 



Xd 



• 1 1 

a i a jXij 



where 



Xij = / Cij(T)dT. 

Jo 



(7) 



(8) 



Here, the 'basic' correlation functions CV,-(t) arc defined 
as: 

Cij(r) = (9) 
\ ((Mr)A 3 {Q)) + (Mr)MO))) - (M0))(M0)) ■ 



Next, we note that multiplication of the correlation 
function Cq{t) by an arbitrary constant factor simply 

corresponds to multiplying the operator O by the square 
of that constant, and we should therefore restrict the dis- 
cussion to normalized correlation functions. This is done 
by the natural requirement that the value of the correla- 
tion function at r = 0, namely Cq(0), be one. In terms 
of the coefficients a,-, this translates to the condition: 



(10) 



where 



Vij = <7«(0) = - [{AiAj) + (AjAi) ) - (Ai){Aj) . (11) 

Note that unlike the matrix elements Xij which depend 
on the long-time behavior of the system, i.e., on the entire 
spectrum of the Hamiltonian, the matrix elements rjij are 
ground-state properties. 

Interestingly, the above normalization condition can be 
expressed as an inner product in this 'space of basic oper- 
ators' {^4j}- For any two arbitrary measurable operators 
A and B, this inner product is simply the equal-time 
covariancc: 

A*B = ±((AB) + (BA))-(A)(B), (12) 
from which the norm 

||i|| = (A*i) 1/2 =((i 2 )-(i) 2 ) 1/2 (13) 

is immediately derived. 

Denoting the vector of the coefficients of O by 
a. = (ai, . . . , Qjj, . . . , cxm), our problem translates to 
maximizing the quantity (a|^|a), where x is th e positive 
definite matrix whose ij-th entries are Xij ; supplemented 
by the normalization condition ||0|| = 1 which translates 
to 



(a|?7|a) = 1 



(14) 



where r\ is the equal-time covariancc matrix whose entries 
are, analogously, rjij. It too is a positive definite matrix. 
Both matrices r\ and % fall under the category of Gra- 
manian matrices for which each entry can be viewed as 
an inner product of two elements chosen from a set of M 
elements. As we shall see later, the two sets of entries, 
Xij and rjij have something else is common: both can be 
very easily measured within QMC simulations. 



B. Projecting out conserved quantities 

Before moving on to the maximization procedure how- 
ever, there is another delicate issue that needs to be ad- 
dressed: It may very well be the case that the Hamilto- 
nian governing the physical system in question has a set 
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of conserved quantities associated with it (one of whom 
would usually be the Hamiltonian itself). Let us denote 
the operators corresponding to these quantities by the set 
{Bk} with k = l..N c where N c is the 'number of linearly 
independent constraints' or equivalently the 'number of 
conserved quantities', and rewrite each of those if possible 
as a linear combination of the set of the basic operators 
[Ai] that comprise 6 [11]: 

B k =J2ti k) A, (15) 

(k) 

where (3\ are real- valued coefficients. The normalized 
correlation functions corresponding to these conserved 
quantities are all simply Cg k (r) = 1, i.e., they are con- 
stant in imaginary time. 

It is important to note then that maximizing \o with- 
out taking these operators into account is guaranteed to 
produce the 'optimized' correlation function Cq(t) = 1, 
with the maximal value of \d = P corresponding to an 
operator which is an arbitrary linear combination of these 
conserved quantities. This of course is a situation that 
one would wish to avoid, as the gap could not be ex- 
tracted from a constant correlation function. It is there- 
fore necessary to 'remove' the above conserved quantities 
from the optimized operator O. Interestingly, in terms 
of the newly defined inner product, this condition may 
be formulated very naturally by requiring that O be or- 
thogonal to, or in other words uncorrelated with, each 
of the operators corresponding to the various conserved 
quantities, namely by requiring that O * Bk = for each 
k. In vector notation, these requirements translate to: 

(oc\r ) \(3^)=0. (16) 

The careful reader will notice that there is a certain 
subtlety associated with the condition Eq. (jTHJ) , which is 
important to address. As discussed earlier, it is crucial 
for the extraction of the gap that our system be strictly in 
its ground state, as we shall assume it is for all practical 
purposes. In the ground state, the conserved quantities 
{Bk} do not fluctuate. They obey 

(Bt)-(B k ) 2 = 0. (17) 

In terms of the equal-time covariancc matrix r] (and 
also in terms of the matrix x) this translates to 
(/3 (fe) |T7|/3 (fe) ) = 0. The above equation implies that r] 
is in practice no longer strictly positive definite but only 
positive semidefinite and the set of vectors {/3^} spans 
its kernel (and also that of x)- I n this case, it would 
be impossible to require that the condition Eq. p^|) be 
satisfied. For the optimized operator O to be orthogonal 
to those, we must restrict O to the subspace orthogonal 
to the kernel of rj, that is, to require that the vector a 
satisfies the amended conditions 

(a|/3«)=0, (18) 

for all k = l..N c . 



In passing, we note the following: It may very well be 
the case that one would not be aware of all the conserved 
quantities associated with a given Hamiltonian and that 
can be given as linear combinations of the basic opera- 
tors. Therefore, in the course of optimizing the corre- 
lation function for the extraction of the gap, the to-be- 
optimized operator O will not be orthogonal to all con- 
served quantities. In this case, the maximization of Xo 
will yield some linear combination of operators from the 
set {Ai} which would again give Cq = 1. The resulting 
operator will produce a correlation function with which 
one could not extract the gap, as the latter would just 
be a constant. The procedure would however yield the 
unknown conserved quantity. Put differently, the maxi- 
mization process described above may be used to detect 
unknown conserved quantities of the system (although 
this would probably require very accurate measurements 
and may therefore turn out to be a numerically very de- 
manding task). Once all the conserved quantities are 
revealed, they may then be removed from O by suitable 
orthogonality conditions. 



C. The optimization process 

At this stage, we can reformulate the problem at hand 
in a purely mathematical form: Given two symmetric 
M x M square matrices x an d V an( i a se t °f vec ~ 
tors (/3 (1) ,...,/3 (fc) ,.--,/3 (Ar<:) ) with N c < M, find the 
vector a that maximizes (a|^;|a) given the constraints 

(a|rj|a) = 1 and (e*|/3 (fe) ) = for all k = l..N c . 

The solution to the above problem is easily obtainable 
and is given by the following prescription: 

• Using a Gram-Schmidt orthonormalization process, 
find a set of (AI — N c ) orthonormal vectors {/3{l } 
with k = 1..M —N c , that span the subspace orthog- 
onal to that spanned by the set {/3^ k '}. Construct 
then the (M — N c ) x M matrix Pj_ whose rows are 

the vectors {f3^}. 

• Define the 'reduced' matrices f) — P±r)P± and 
X = P±xPl These are just rj and x with the 
subspace spanned by the set {/3^ fc '} removed. In 
our case, this ensures that fj is a positive definite 
matrix. The dimensionality of each of the new ma- 
trices is (AI - N c ) x (AI - N c ). 

• Now the problem reduces to finding a vector a that 
maximizes: 



The above expression would be maximal if we set 
|a) to be the eigenvector of J7 -1 / 2 xrj -1 / 2 belonging 
to the largest eigenvalue. 
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• Switching back to the full Hilbcrt space, the vector 
\a) =P[f)-^ 2 \a) (20) 
is then the solution to our problem. 

D. Practical guidelines and cost of the procedure 

Taking time-correlation measurements of the opti- 
mized composite operator O = Yli=i a iAi using the op- 
timal set requires first knowing the values of these 
coefficients up to an acceptable statistical error. In the 
previous section, we saw that the values of these coeffi- 
cients are given as functions of the matrix elements \ij 
and rjij and these need to be determined beforehand. 
Thus, measuring the imaginary-time correlations of O 
requires that the QMC simulation has two 'phases' as far 
the gap calculations are concerned: In the first phase, the 
optimal set {on\ needs to be determined by performing 
measurements of Xij and f]ij until the desired accuracy is 
reached. At the end of this phase, the set of parameters 
{at} is calculated according to the prescription given in 
Sec. IIII CI Thus, during the first phase of the simulation 
no imaginary-time correlations are measured; the actual 
measurements of the correlation function corresponding 
to the optimal operator are performed in a second phase 
of the simulation. 

Measurements of M 2 physical quantities during the 
first phase of the simulation may seem a bit costly at first 
due to the fact that the number of independent 'basic' op- 
erators in a given problem usually scales like the number 
of particles in the system TV, that is, M 2 ~ TV 2 . However, 
we note that this seemingly high cost is compensated by 
the fact that the number of operations needed for this 
measurement process does not scale with, and in fact is 
independent of, the inverse temperature /3. For ground- 
state measurements, the appropriate inverse-temperature 
P normally grows polynomially or even exponentially 
with TV. This is because of the condition of f3AEi 3> I 
which needs to be maintained while the gap to the first 
excited state AEi usually decreases at least polynomi- 
ally fast in TV. It is therefore plausible to assume that 
TV/ 2 < TV/3 and so, the procedure of calculating the above 
matrix entries comes at a rather low price: It requires 
less than 0(N/3) operations. 

In practice, finding the kernel (i.e., the subspace 
spanned by the vectors representing conserved quanti- 
ties) of the equal-time covariance matrix 77 numerically 
may turn out to be a rather difficult task especially for 
large system sizes. This is because of the statistical errors 
associated with the measured matrix elements of r\ which 
may eventually lead to negative eigenvalues, despite the 
fact that 77 should be strictly positive definite. The ex- 
istence of negative eigenvalues implies that the errors 
and corresponding negative eigenvalues are comparable 
in size. Existence of negative eigenvalues also has harsh 
consequences as far as the maximization procedure de- 
tailed in the previous section is concerned, as the process 



requires that r\ be a positive-definite matrix. Therefore, 
a practical resolution of the above difficulties would be to 
simply add the cigenstatcs associated with the negative 
eigenvalues to the set of vectors {(3^} thereby ensuring 
that the subspace spanned by the remaining eigenstates 
of rj has strictly positive eigenvalues as required by the 
maximization process. As we shall later see, this solution 
works very well in practice. 

Moving on to the second phase of the simulation, the 
operator O is constructed using the optimized set of pa- 
rameters whose values were set at the end of the first 
phase of the simulation, and its time correlations are 
measured. The construction of the composite operator O 
consists of 0(Nj3) operations corresponding to its evalu- 
ation (M ~ TV operations) at each "time-slice" of which 
there are of the order of (3. Since the time correlations of 
only one operator are being measured, calculation of the 
actual correlation data requires of the order of 2/3 log j3 
operations if one uses the Fast Fourier Transform algo- 
rithm to compute them. This procedure may therefore 
be considered cheap computationally as well. 

Since the measurement of the optimized correlation 
function requires two phases of the simulation, it may 
seem that it also demands more computation time as 
compared to correlation function measurements that do 
not require optimization. In practice however, one finds 
that this is not the case. When using an optimized cor- 
relation function, two effects come into play. First, the 
slowest-decaying exponent of the correlation function be- 
comes dominant at much shorter imaginary times and 
second, the correlation function being 'lifted' above the 
horizontal axis results in a much less noisy curve. In 
practice, as we shall see in the next section, these effects 
translate to an overall effect of a much shorter needed 
computation time. 

E. The two-operator case 

The simplest nontrivial example for a somewhat- 
optimized measurable operator one can think of and 
which can be given a closed-form expression, is one in 
which the operator to be optimized is a combination of 
only two other operators, namely, O = a\0\ + 
where 0\ and O2 correspond to two measurable quan- 
tities, each corresponding to some linear combination of 
the basic operators. Assuming also that the Hamilto- 
nian is a linear combination of the two operators, namely 
H = (3\0\ + P2O2, the condition (a\f3) = 0, along with 
the normalization condition of O and the fact that the 
vector (/3i,/32) spans the kernel the associated 2x2 ma- 
trix T7, yields the answer: 

( , / -/3i/3 2 (-/3 2 ,/3 1 ) 

(ai ' a2 ^V M+M ' { ) 

where no maximization of (a|x|ci) is needed. In this 
case, the simulation need not be split into two phases in 
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order to determine the optimal coefficients. 

IV. ILLUSTRATIVE EXAMPLES 

In what follows, we demonstrate the effectiveness of 
using an optimized imaginary-time correlation function 
to extract the excitation gap. We do this by considering 
several problems taken from the field of Quantum Adi- 
abatic Computation, in the context of which the Quan- 
tum Adiabatic Algorithm (QAA) [2l| has been devised 
to solve hard optimization problems efficiently on a quan- 
tum computer. 

Within the framework of this approach, the efficiency, 
or complexity, of the QAA for a given input problem 
is often studied by analyzing the behavior of the exci- 
tation gap of a one-parametric family of Hamiltonians 
that forms a linear interpolation between an easily solv- 
able transverse-field Hamiltonian Hd (commonly referred 
to as a 'driver' Hamiltonian) and a diagonal 'problem' 
Hamiltonian H p whose ground state encodes the solution 
to the optimization problem. Put explicitly, the linear in- 
terpolation is 

H(s) = sH p + {l-s)H d , (22) 

where s G [0, 1]. In these problems, the gap usually needs 
to be calculated for several values of s where the objective 
is to find the minimal gap among these (the reader is 
referred to Refs. pTI42^ | for a more detailed description 
of the process). 

Here, we shall illustrate the advantages that come 
with using the optimization method described in previ- 
ous sections by considering several typical instances of 
a specific optimization problem of the "constraint sat- 
isfaction" type known as l-in-3SAT (for a description 
of the problem see, e.g., Refs. [22| and [IBj]), in which 
the Hamiltonian is a sum of L three-local Hamiltonians, 
Hp = ^2^ = i H a , where each term in the sum involves 
three spins picked randomly from a pool of N spins. Each 
local Hamiltonian H a is given in this problem by the ex- 
pression: 

H a = -(5-cr* -a z „ -o z n (23) 

1 „z „Z 1 Z Z 1 z z 1 o„z z z \ 

+ a a 1 G a, + °a 2 a a 3 + G3 o 1 + 60 ai ° a 2 ° a 3 J > 

where a; for i = 1,2,3 label the participating spins and 
erf is the z-component Pauli matrix acting on spin i. The 
second part of the Hamiltonian is the driver Hamiltonian 
Hd- It is a simple transverse-field Hamiltonian and is 
given by: 

^=-7jE<> ( 24 ) 

i=i 

where erf is the x-component Pauli matrix acting on spin 
i. 



A. A 16-spin system 

Here, we analyze an instance of the l-in-3SAT problem 
in which the number of spins in the system is N = 1 6 
and the number of clauses, each involving a randomly- 
chosen triplet of spins, is L = 13. The relatively small 
size of the system will allow us to compare the QMC- 
based extracted gap against the corresponding exact- 
diagonalization result. The chosen inverse temperature 
is j3 = 128 (which, as we shall see obeys j3AEi 3> 1) and 
the value chosen for the parameter s in this example is 
s = 0.5, which puts it very close to the location of the 
minimum gap. 

To illustrate the effectiveness of using an optimized 
correlation function, we will also compare our results 
against those obtained from a 'partially-optimized' cor- 
relation function based on an optimization with respect 
to only two coefficients, and a 'non-optimized' correla- 
tion function whose coefficients are not optimized by any 
means but are chosen randomly instead. 

The specific QMC method we use here to measure the 
excitation gap is known as the stochastic series expansion 
(SSE) algorithm (HH^. This method involves a Taylor 
series expansion of the partition function Z = Trfe - ^] 
and uses a discrete representation of continuous imag- 
inary time. Similarly to current path integral formu- 
lations in continuous imaginary time |28l430l |. this dis- 
cretization does not introduce errors into the algorithm. 

Within the SSE scheme applied to the problem at 
hand, the 'natural' operators to measure are the -/V non- 
diagonal operators Ai = — 1(1 — s)erf with i = 1..N and 

the L diagonal operators Ai = sH a with a = 1..L where 
i = N + a. In the current example this amounts to a 
total of M = N + L = 29 basic operators. The to-be- 
optimized operator O will thus be a linear combination 
of those. Associated with the model studied here is only 
one conserved quantity - the energy of the system (hence 
the number of conserved quantities is N c = 1). It corre- 
sponds to the set of coefficients /8- = I (with i = 1..M). 

To obtain the optimal operator O, we calculated in 
the course of the simulation the matrix entries Xij an d 
rjij each corresponding to an easily measured quantity 
[see Eqs. © and (fTTj) ]. The expressions for the various 
matrix elements as they are implemented within the SSE 
algorithm are derived in Appendix [X] for the convenience 
of the reader. 

As a next step, we followed the procedure outlined in 
Sec. IIIICI and diagonalized the matrix if in order to find 
its kernel. As expected from positive semi-definite ma- 
trices, the resulting eigenvalues were all found to be pos- 
itive. One particularly small eigenvalue was also found 
and was immediately identified as the 'zero' eigenvalue 
corresponding to the one conserved quantity - the energy. 
Numerically, the eigenvalue was found to be of the order 
of I0~ 7 , whereas for comparison, the next smallest eigen- 
value turned out to be about 10 3 times greater. We also 
found as expected, that both matrices shared the same 
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10 15 20 
operator index 



FIG. 2: (Color online) Optimized (red triangles), 'partially- 
optimized' (green squares) and 'non-optimized' (i.e., a 
randomly-generated operator - blue circles) coefficients val- 
ues of the operator O whose different-time correlations are 
used to extract the gap. The values shown are normalized 
according to ||0|| = 1 or equivalently {a\rj\a) — 1. The hori- 
zontal line indicates the indices of the basic operators as they 
are defined in the text. 

kernel corresponding the eigenvector (1, 1, . . . , 1) which 
represents the Hamiltonian of the system. After pro- 
jecting out the kernel of r\ from the two matrices (i.e., 
reducing the matrices to f) and x), we calculated the 
coefficients {a{\ by maximizing the expression given in 
Eq. (fT9|) to obtain the optimal operator O. These coeffi- 
cients are plotted in Fig. [2] (marked by triangles). As the 
figure indicates, the various coefficients on turn out to be 
quite different from one another, revealing the nontriv- 
ial complexity of the optimal operator. The correlation 
function itself is plotted in Fig. [3] accompanied by a linear 
fit with which the excitation gap is extracted. 

For comparison, we have also considered a much sim- 
pler set of operators, namely 0\ = (1 — s)Hd = J2iLi 
and Oi = sH p = $3»=jv+i Ai> to construct the operator 
O from. For this choice of operators the coefficients of 
the Hamiltonian correspond to = 13% = 1. Plugging 
these into Eq. pTj) . we end up with ot\ = — ct2 and the 
(unnormalized) 'partially-optimized' operator is: 



O p o . = sH p - (1 - s)H d . 



(25) 



Measurement of the correlation function corresponding 
to this operator is also depicted in Fig. [3] (square data 
points), and the corresponding coefficients arc plotted in 
Fig. [2] for comparison. 

Figs. [2] and [3] also show the (normalized) coefficients 
and correlation function of a 'non-optimized' operator 
constructed by random assignments of the various coef- 
ficients taken from a uniform distribution in the range 
[—1,1] (these are marked by circles in the two figures). 

A sidc-by-side comparison of the three correlation 
functions plotted in Fig. [3J clearly indicates that while 
the non-optimized correlation function is very noisy and 



Operator-type 


Gap value and its error 

(the relative error is 
shown in parentheses) 


Deviation from 
ED (in standard 
deviations) 


fully-optimized 
partially-optimized 
non-optimized 


0.0747 ± 0.0006 (0.8%) 
0.076 ± 0.005 (6.6%) 
0.053 ± 0.02 (38%) 


0.38 
0.31 
1.07 



TABLE I: Calculated excitation gaps for the 16-spin sys- 
tem as they were extracted from the three tested correlation 
functions, and their deviation from the exact-diagonalization 
(ED) value measured in number of standard deviations. As 
the table indicates, the optimized correlation function yields a 
much better prediction of the gap than the partially- and non- 
optimized correlation functions. The exact diagonalization 
value is AE\ — 0.07447 and the chosen inverse-temperature 
here is /3 = 128. The requirement that the system be in its 
ground state is therefore satisfied as (3AE « 10. The errors 
reported here were determined by standard least-squares fit- 
ting analysis on the correlation-function data. 

simply does not allow for any serious calculation of the 
gap at least for this chosen running time of the simula- 
tion, the two other optimized correlation functions are 
much less noisy: both exhibit a straight-line behavior at 
some point in imaginary time and enable the extraction 
of the gap. That being said, it is evident from the fig- 
ure that the dominance of the slowest-decaying exponent 
is apparent in the fully-optimized correlation function at 
much shorter times than it is for the partially-optimized 
correlation function. Moreover, the fully-optimized cor- 
relation function also has a superior signal-to-noise ratio 
at all times. 

Calculation of the excitation gap by the linear 
fits to the logarithm of the three (fully-optimized, 
partially-optimized and non-optimized) correlation func- 
tions yields the results summarized in Table |U As 
the table indicates, for the parameters chosen here, 
a linear fit of the fully-optimized correlation function 
clearly yields the most accurate result among the three, 
with a very small error (of 0.0006) and the exact- 
diagonalization value falling well within the error bar 
of the extracted value (0.38 standard deviations). The 
partially-optimized correlation function yields slightly 
poorer results (an error of 0.005 and with the exact value 
being 0.31 standard deviations away), although one could 
argue that such 'partial' optimization might be enough 
for certain purposes. In this case, with essentially only 
one adjustable parameter, the partially-optimized cor- 
relation function is only slightly worse than the fully- 
optimized one. The non-optimized data on the other 
hand produced much larger errors (0.02) and a rather 
poor fit, as one would expect from such noisy correlation 
function. 



B. Larger spin systems 

In what follows, we present the results of an analy- 
sis similar to the one performed in the previous section 
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System 
size 


Operator-type 


Gap value 
and its error 


Relative 
error 


TV = 48 


fully-optimized 
partially-optimized 


0.065 ± 0.0015 
0.064 ± 0.005 


2.3% 
7.8% 


TV = 64 


fully-optimized 
partially-optimized 


0.051 ±0.001 
0.052 ± 0.007 


2% 
13% 



FIG. 3: (Color online) Optimized (red triangles), 'partially- 
optimized' (green squares) and 'non-optimized' (blue circles) 
correlation functions as functions of imaginary time on a 
log-linear scale for an instance of a l-in-3SAT problem with 
TV = 16 interacting spin- 1/2 particles. As the figure indicates, 
the optimization process substantially reduces the contribu- 
tions from higher-order coefficients leaving even in relatively 
short times only the slowest-decaying exponent. The resulting 
optimal correlation function was easily fit with a linear func- 
tion at the appropriate region producing an accurate predic- 
tion of the excitation gap (the various linear fits are the solid 
lines). The partially-optimized correlation function is also a 
vast improvement over the non-optimized correlation function 
which in turn yields results of a much poorer quality. For the 
latter correlation function, there seems to be no region where 
a linear behavior is obvious, as higher-order contributions are 
evident throughout the (also very noisy) examined region. 

but applied to larger spin systems. For these, the ex- 
pected gap is much smaller and is also naturally more 
difficult to extract. Here we shall consider two random 
instances of the l-in-3S AT type corresponding to systems 
with TV = 48 and TV = 64 spins, and with L = 38 and 
L = 51 clauses, respectively. While in these examples a 
comparison with exact-diagonalization results is unavail- 
able, it is nonetheless advantageous to make a compari- 
son between the fully-optimized, partially-optimized and 
non-optimized correlation functions. 

The correlation functions obtained for the TV = 48 sys- 
tem (the adiabatic parameter s and inverse temperature 
j3 were chosen to be s = 0.5 and (3 = 256) are shown in 
Fig. @] along with corresponding linear fits from which the 
gaps are then extracted. The optimized correlation func- 
tion was obtained via the maximization of the integrated 
susceptibility with respect to all of the M = 48 + 38 = 86 
coefficients on. In this specific example, the subspace 
spanned by five eigenvectors corresponding to five nega- 
tive eigenvalues that were found at the end of the first 
phase of the simulation have been discarded. The sec- 
ond correlation function is the partially-optimized one 
corresponding to the operator O p o . given in Eq. (|2~5T) . 
and the third correlation function tested here is the non- 
optimized one, based on randomly-assigned coefficients. 

As in the previous example, it is evident from Fig. 0] 
that the fully-optimized correlation function is much 



TABLE II: Calculated excitation gaps for the 48- and 64- 
spin systems as they were extracted from the three tested 
correlation functions and their relative errors. As the table 
indicates, the optimized correlation function yields a much 
more accurate prediction of the gap than the partially- and 
non-optimized correlation functions for both sizes. 

less noisy and therefore produces a much more accu- 
rate reading of the gap with the smallest uncertainty 
AEi = 0.065 ±0.0015 (2.3% relative error). For compar- 
ison, the partially-optimized correlation function yielded 
AE 1 = 0.064 ± 0.005 (7.8% relative error). The non- 
optimized correlation functions turned out to have very 
large statistical errors and is therefore completely unre- 
liable for the extraction of the gap. These results are 
summarized in Tabic [II] 

Similarly to the TV = 48 case analyzed above, we 
present in Fig. [5] the correlation functions obtained in 
the analysis of a 64-spin system with 51 clauses and a 
total of M = 115 real coefficients (here, s = 0.52 and 
j3 = 256). In this example, a 23-dimcnsional subspace has 
been discarded in the maximization procedure, leaving 
a 92-dimensional parameter-space with respect to which 
maximization is performed. As in the previous exam- 
ples, here too the non-optimized, i.e., random, correla- 
tion function is practically useless for obtaining the gap, 
whereas the partially-optimized and fully-optimized cor- 
relation functions are much more suitable for the task. 

As expected, the fully-optimized correlation function 
is significantly less noisy than the partially-optimized 
one giving a much more accurate value for the gap: 
AE ± = 0.051 ± 0.001 (2% relative error) for the fully- 
optimized versus AEi = 0.052 ± 0.007 (13% relative 
error) for the partially-optimized one. Again, the non- 
optimized correlation function is much too noisy for any 
reliable reading of the gap (see also Table Hlj) . 



V. SUMMARY AND CONCLUSIONS 

In this paper, we demonstrated the effectiveness of 
calculating and utilizing optimized imaginary-time cor- 
relation functions toward the extraction of excited-state 
information within quantum Monte Carlo simulations. 
We have given a prescription for optimizing the operator 
whose imaginary-time correlation function would be best 
suited for gap calculations. The optimization is based 
on maximizing the integrated susceptibility of the op- 
erator whose time-correlations are to be measured. Wc 
have also illustrated the benefits associated with evalu- 
ating optimized correlation function compared to other 
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FIG. 4: (Color online) Optimized (red triangles), 'partially- 
optimized' (green squares) and 'non-optimized' (blue circles) 
correlation functions as functions of imaginary time on a log- 
linear scale for an instance of the l-in-3SAT problem with 
N = 48 interacting spin-1/2 particles. The figure shows that 
the optimal correlation function is easily fit with a linear func- 
tion at the appropriate region producing an accurate predic- 
tion of the excitation gap (the various linear fits are the solid 
lines). The partially-optimized correlation function is also a 
vast improvement over the non-optimized correlation function 
which in turn yields results of a much poorer quality. 



was presented here may of course be subjected to other 
more sophisticated methods aimed at obtaining a fuller 
picture of the spectrum of the Hamiltonian of the system, 
i.e., energies of more excited states (see e.g., Ref. [iH - 
[lH). These methods of analysis can readily be applied 
to the correlation function obtained in the process de- 
scribed here. 

The results of the analysis performed on the exemplary 
problems presented in the Sec. IIVI show that optimiza- 
tion of the correlation function results in both smaller 
statistical errors of the correlation function and also a 
dominance of the slowest-decaying exponent at shorter 
imaginary times. This implies that in order to obtain de- 
cent results for the gap from non-optimized correlation 
functions, a substantially longer running time of the sim- 
ulation is needed (and perhaps also a larger value of the 
inverse-temperature /?). The results therefore establish 
the importance of applying at least some optimization to 
measured correlation functions, as this may have consid- 
erable effects on the ability to extract the gap and on the 
accuracy of the obtained value. We therefore expect that 
the procedure presented in this manuscript will be useful 
wherever gap calculations are needed numerically. 
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Appendix A: Stochastic Series Expansion (SSE) 
measurements of the matrices rj and x 



FIG. 5: (Color online) Same as Fig. [4] but for a system with 
TV = 64 spins. Here too the optimal correlation function 
(red triangles) is much less noisy than the other partially- 
optimized (green squares) and non-optimized (blue circles) 
correlation functions, eventually leading to a more accurate 
prediction of the system gap. Of the two latter correlation 
functions, the non-optimized one is particularly noisy in this 
example and yields an extremely poor estimation of the gap. 

non-optimized functions, and confirmed numerically that 
determining the gap to the first excited state from opti- 
mized functions is considerably more accurate than cor- 
responding results obtained from non-optimized correla- 
tion functions. Wc have also commented on the relatively 
low cost of the procedure. 

While in this study the main focus was on the extrac- 
tion of the gap to the first excited state, it should be 
noted that the optimization procedure whose derivation 



In what follows we derive the explicit expressions 
needed for the measurements of the matrix entries rjij 
and Xij within the framework of the stochastic series ex- 
pansion (SSE) algorithm. As discussed in Sec. IIII1 cal- 
culation of the optimal correlation function requires the 
evaluation of the set of coefficients {o^} which is obtained 
by an algebraic manipulation of the matrix elements \ij 
and rj^ as they are defined in the main text [see Eqs. ((SJ) 
and (fTT]) ]. These are expressed in terms of the 'basic' 
operators of the QMC technique being used. Within the 
SSE, the basic operators that are most easily measured 
are the so-called 'bond' operators which we denote here 
by Ai. These are the local operators that comprise the 
Hamiltonian, namely: 

M 

H = J2 A i> ( A1 ) 

i=i 
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where each of the bond operators Ai acts only on a lim- 
ited number of particles. The reader is referred to Sec. lIVI 
for the bond operators of the l-in-3SAT problem. Within 
the SSE, the expectation values of the bond operators 
are obtained by simply counting their occurrences in the 
'operator sequence' that defines the instantaneous con- 
figuration of the system (for a more detailed discussion 
of the SSE technique, see, e.g., HE H3): namely: 

(A i ) = -~(N i ), (A2) 

where Ni is the number of times the operator Ai appears 
in the sequence. Similarly, it is easy to show that 

r fi i 

j &T{A tl { T )A l2 {0)) = piiNiM - Si^iNi,)) (A3) 

where i\ and ii are two arbitrary bond indices [26]]. From 
the above equation, it follows that the matrix elements 



of x are simply given by: 

Xij = ^ m,N i2 ) - Si^iNiJ - (N n )(N l2 )) . (A4) 

In addition, expectation values of products of bond op- 
erators are obtained using: 

(H il H i2 ) = jp((n-l)N il2 ), (A5) 

where Ni 12 denotes the number of ordered subsequences 
Ai x Ai 2 in the operator sequence and n = X)f=i Ni is the 
total number of operators in the sequence. The matrix 
elements of r\ are therefore similarly given by the expres- 
sion: 

Vij = i (\{{n - l)(N ll2 + N i21 )) - (N^nJ) .(A6) 
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For the discussion here we shall consider the canonical 
ensemble scheme, although the following may just as well 
be applied to the grand-canonical ensemble. 
Of course, in some cases where critical slowing down is 
unavoidable, especially in the vicinity of first order phase 
transitions, relaxation to the ground state may take an 
exponentially-long amount of time. 

Obviously, one could choose to work with only a subset 
of the available operators or rather with combinations of 
them, thereby reducing the total number of operators M. 
It should be noted that in principle there could be a situa- 
tion where a conserved quantity could not be constructed 
from the set of basic operators. 



